Optimization hardness in biological transients¶

Preamble¶

  • After opening this notebook, run the cells below before anything else. These will import libraries and define functions.

  • This notebook requires the following packages, which are pre-installed on Google Colab.

    • Python 3
    • NumPy
    • SciPy
    • Matplotlib
In [ ]:
## Preamble / required packages
import numpy as np
from IPython.display import Image
np.random.seed(0)

## Import local plotting functions and in-notebook display functions
import matplotlib.pyplot as plt
from IPython.display import Image, display
%matplotlib inline
## default first plot color is black
plt.rcParams['axes.prop_cycle'] = plt.cycler(color=['k'])  # Set default plot color to black
plt.rcParams['axes.prop_cycle'] = plt.cycler(color=plt.cm.tab10.colors)  # Reset to default colors
plt.rcParams['legend.frameon'] = False

## Disable warnings
import warnings
warnings.filterwarnings('ignore')

# !apt-get install nvidia-cuda-toolkit
# !pip3 install numba
## These lines allow Google Colab to import and work with external files.
# %rm -rf illotka
# !git clone https://github.com/williamgilpin/illotka.git
# %cd illotka
# import importlib.util, sys
# module_path = "/content/illotka/base.py"
# spec = importlib.util.spec_from_file_location("base", module_path)
# base = importlib.util.module_from_spec(spec)
# sys.modules["base"] = base
# spec.loader.exec_module(base)

# import requests
# url = "https://raw.githubusercontent.com/williamgilpin/illotka/main/base.py"
# resp = requests.get(url)
# resp.raise_for_status()
# with open("base.py", "w") as f:
#     f.write(resp.text)

# url = "https://raw.githubusercontent.com/williamgilpin/illotka/main/utils.py"
# resp = requests.get(url)
# resp.raise_for_status()
# with open("utils.py", "w") as f:
#     f.write(resp.text)


from base import *
from utils import *

# Autoreload
%load_ext autoreload
%autoreload 2
print("Loading complete")
Loading complete
In [ ]:
 

Revisiting stability vs. complexity¶

  • Biological systems,

  • May 1972:

In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 

The generalized Lotka-Volterra model¶

We consider random ecosystems given by the generalized Lotka-Volterra equation,

$$ \frac{dN_i}{dt} = N_i \left( r_i + \sum_{j=1}^N A_{ij} N_j \right) $$

where $N_i$ is the population of species $i$, $r_i$ is the intrinsic growth rate of species $i$, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$.

To study the behavior of this model, we consider the case where $r_i \sim \mathcal{N}(0,1)$, and the matrix $A$ has the form

$$ A = Q - d\, I $$ where $I \in \mathbb{R}^{N \times N}$ is the identity matrix, $Q_{ij} \sim \mathcal{N}(0,1)$ and $d$ is a constant density-limitation.

In [2]:
# Set the simulation parameters
tmax = 20000 # Maximum integration duration
n = 200 # Number of species
d = 12.5 # Density-limitation

## Initialize the model
eq = GaussianLotkaVolterra(n)
eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
eq.r = np.random.normal(size=n)

## Initial conditions
ic = 0.4 *np.random.uniform(size=n)

## Numerically integrate the model
t, y = eq.integrate(tmax, ic)

plt.figure(figsize=(7, 4))
plt.semilogx(t, y, color="k", lw=1.5, alpha=0.3);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, np.max(t))
plt.ylim(0, None)
Out[2]:
(0.0, 0.6341573011172086)
No description has been provided for this image

How many species survive?¶

In [3]:
nonzero_remaining = np.sum(y[-1] > 1e-12) # Precision floor
print(f"{nonzero_remaining} species out of {n} remain")
106 species out of 200 remain

We repeat the experiment many times¶

In [5]:
n_replicates = 20
tmax = 1000

number_survive = []
for i in prange(n_replicates):
    progress_bar(i, n_replicates) # Print progress bar for loop

    eq = GaussianLotkaVolterra(n)
    eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
    # eq.A = eq.A * np.identity(n) # Wipe out all off-diagonal elements
    eq.r = np.random.normal(size=n)
    t, y = eq.integrate(tmax, ic)

    n_survive = np.sum(y[-1] > 1e-12)
    number_survive.append(n_survive)

plt.hist(number_survive)
plt.xlabel("Number of surviving species");
plt.ylabel("Frequency");
[################### ] 

No description has been provided for this image
In [31]:
print(np.mean(number_survive))
print(np.std(number_survive))
103.66
6.259744403727679
In [ ]:
 
In [32]:
print(0.5 * n)
print(np.sqrt(n * 0.5 * (1 - 0.5)))
100.0
7.0710678118654755
In [ ]:
 
In [ ]:
 
In [ ]:
 

Properties of the Random Lotka-Volterra model¶

Proven in Serván et al. 2018.

  1. On average, half of all species survive at steady-state. Across replicates, the number of surviving species follows a binomial distribution with mean $n/2$.
  2. For a given matrix $A$, the steady-state solution is unique (global attractor).
  3. For a given matrix $A$, there always exists a density-limitation $d$ that that leads to a single stable steady-state.

Main assumption: Interactions and growth rates are drawn from a symmetric distribution with finite variance.













How do extreme values affect the dynamics?¶

The generalized normal distribution is parametrized by shape parameter $\beta$, which controls the tail thickness. It is given by

$$ p(x) = \frac{\beta}{2\sigma\Gamma(1/\beta)} e^{-\left(\frac{|x-\mu|}{\sigma}\right)^\beta} $$

  • when $\beta < 1$, the distribution has fatter tails than the normal distribution.
  • when $\beta > 1$, the distribution has lighter tails than the normal distribution.
In [30]:
# Set the simulation parameters
tmax = 10000 # Maximum integration duration
n = 300 # Number of species
d = 18 # Density-limitation
seed = 0 # Seed for random number generator
# seed = None

## Initialize the model
np.random.seed(seed)
eq = GaussianLotkaVolterra(n, random_state=seed, tolerance=1e-10)

## Sample the interaction matrix from a normal distribution
# eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
# eq.r = np.random.normal(size=n)

## Sample the interaction matrix from a modified normal distribution
eq.A = 1.5 * normal_generalized(size=(n, n)) - d * np.identity(n)
eq.r = normal_generalized(size=n)

## Random initial conditions
ic = 0.5 * np.random.uniform(size=n)

## Numerically integrate the model
t, y = eq.integrate(tmax, ic)

## Plot the results
plt.figure(figsize=(7, 4))
plt.semilogx(t, y, color="k", lw=1.5, alpha=0.2);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, 1000)
plt.ylim(0, 0.1)


nonzero_remaining = np.sum(y[-1] > 1e-12) # Precision floor
print(f"{nonzero_remaining} species out of {n} remain")
156 species out of 300 remain
No description has been provided for this image
In [ ]:
 
In [ ]:
 

How long does it take to reach the steady-state?¶

Given a steady-state solution $\mathbf{N}^*$, the settling time is the time it takes for the solution to reach a given tolerance $\epsilon$ of the steady-state.

$$ \argmin_{t} \bigg( \dfrac{\| \mathbf{N}(t) - \mathbf{N}^* \|}{\| \mathbf{N}(t) \|} < \epsilon \bigg) $$

In [41]:
def compute_settling_time(traj, tvals=None, eps=1e-6):
    """
    Compute the settling time from a multivariate time series. 
    If the trajectory does not reach the steady-state within the integration time, 
    return a NaN.
    """
    final_value = traj[-1]
    scaled_distance = (np.linalg.norm(traj - final_value, axis=1)) / np.linalg.norm(traj, axis=1)
    final_index = np.where(scaled_distance < eps)[0]
    if len(final_index) == 0:
        return np.nan
    else:
        return final_index[0]



## Initialize the model
np.random.seed(seed)
eq = GaussianLotkaVolterra(n, random_state=seed, tolerance=1e-10)

## Sample the interaction matrix from a modified normal distribution
eq.A = 1.5 * normal_generalized(size=(n, n), beta=1) - d * np.identity(n)
eq.r = normal_generalized(size=n)

## Random initial conditions
ic = 0.5 * np.random.uniform(size=n)

## Numerically integrate the model
t, y = eq.integrate(tmax, ic)

tind = compute_settling_time(y)
print(f"Settling time: {t[tind]:.2f}")
nonzero_remaining = np.sum(y[-1] > 1e-12) # Precision floor
print(f"{nonzero_remaining} species out of {n} remain")
Settling time: 585.71
146 species out of 300 remain
In [ ]:
 

Hierarchical structure in the Lotka-Volterra model¶

Trophic interactions can produce extremal values in the Lotka-Volterra model. We consider sampling a family of matrices of the form

$$ A = P^\top (A^{(0)} - d\, I) P + \epsilon \, E, $$

  • As before,$A^{(0)}_{ij}, r_i \sim \mathcal{N}(0, 1)$

  • $P \in \mathbb{R}^{N \times N}$ is a low-rank assigment matrix describing functional redundancy or trophic levels among species.

  • $E \in \mathbb{R}^{N \times N}$ is a matrix of random noise that prevents exact redundancy.

  • $\epsilon \ll 1$ is a small constant that controls the degree to which redundant species deviate from exact redundancy.

In [ ]:
## Loads a figure from the resources folder
Image("./resources/fig_model.jpg", width=1200)
Out[ ]:
No description has been provided for this image

For an analysis of equilibrium properties of a Lotka-Volterra model with hierarchical structure like this, see Tikhonov (2017).

In [ ]:
 
In [ ]:
 
In [ ]:
m_trophs = 3
n_species = 6
tmax = 1e8

interactions_trophic = np.array([
    [0, 0.5, 0],
    [-0.5, 0, 0.3],
    [0, -0.8, 0],
])*50

p_mat = np.array(
    [[1, 0, 0, 0, 0, 0],
     [0, 1, 1, 0, 0, 0],
     [0, 0, 0, 1, 1, 1]]
)


eq = RandomLotkaVolterra(n_species)
eq.r = np.array([0.8, 0.9, 1]) @ p_mat
eq.A = p_mat.T @ (interactions_trophic - 5*np.identity(m_trophs)) @ p_mat #+ 1e-6 * np.random.normal(size=(n_species, n_species))

ic = np.ones(n_species) * 0.1
t, y = eq.integrate(tmax, ic)

plt.semilogx(t, y);
plt.legend(np.arange(n_species), frameon=False);
plt.xlim(1e-3, None);
plt.xlabel("Time");
plt.ylabel("Species abundance");
No description has been provided for this image

Now, let's scale up¶

  • Instead of imposing a specific functional redundancy, we'll sample the assignment matrix $P$ from a family of random low-rank matrices.

  • The small parameter $\epsilon$ controls the degree to which redundant species deviate from exact redundancy.

In [9]:
tmax = 1e8
n_val = 200

ic = np.random.uniform(size=n_val)
eq = RandomLotkaVolterra(n_species=200, eps=1e-3)
t, y = eq.integrate(tmax, ic)
plt.semilogx(t, y);
plt.xlim(1e-3, None)
Out[9]:
(0.001, np.float64(145205.60974303592))
No description has been provided for this image
In [159]:
tmax = 1e8

n_val = 200
ic = np.random.uniform(size=n_val)
settling_times = []
eps_vals = np.logspace(-5, 0, 10)
for eps in eps_vals:
    print(f"eps: {eps}")
    eq = RandomLotkaVolterra(n_val, random_state=0, eps=eps,  sigma=2.0, d=4.5, kfrac=0.2, connectivity=0.05)
    t, y = eq.integrate(tmax, ic)
    tind = compute_settling_time(y)
    settling_times.append(t[tind])


plt.loglog(eps_vals, settling_times, '.')
plt.xlabel("Distance to Exact Functional Redunancy (epsilon)")
plt.ylabel("Settling Time")
eps: 1e-05
eps: 3.5938136638046256e-05
eps: 0.0001291549665014884
eps: 0.0004641588833612782
eps: 0.0016681005372000592
eps: 0.005994842503189409
eps: 0.021544346900318846
eps: 0.07742636826811278
eps: 0.2782559402207126
eps: 1.0
Out[159]:
Text(0, 0.5, 'Settling Time')
No description has been provided for this image

The Lotka-Volterra model as analogue linear regression¶

We consider random ecosystems given by the generalized Lotka-Volterra equation,

$$ \frac{dN_i}{dt} = N_i \left( r_i + \sum_{j=1}^N A_{ij} N_j \right) $$

where $N_i$ is the population of species $i$, $r_i$ is the intrinsic growth rate of species $i$, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$. The steady-state solutions of this equation has the form

  • Without the outer term, the Lotka-Volterra model is a linear time-invariant system.

  • The solution to the Lotka-Volterra model is the solution to the linear system $A \mathbf{N}^* = -\mathbf{r}$.

  • The outer term adds a positivity constraint. In practice, this leads to sparsity (exclusion) in the steady-state solution $\mathbf{N}^*$.

In [ ]:
eq = RandomLotkaVolterra(n_val, random_state=3, eps=1e0, sigma=0.5, d=4.5, kfrac=0.0, connectivity=0.05)
ic = np.random.uniform(size=n_val)
t, y = eq.integrate(tmax, ic, tol=1e-10)
plt.plot(y);
plt.xlim(1e-3, None)
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.show()
No description has been provided for this image
In [162]:
## Perform unconstrained linear regression
plt.figure(figsize=(6, 3))
plt.plot(np.linalg.solve(eq.A, -eq.r), y[-1], '.')
plt.xlabel("Linear regression solution")
plt.ylabel("Steady-state of Lotka Volterra model")
plt.title("Steady-state solution vs. final time-series value")
plt.show()

## Perform constrained linear regression
from scipy.optimize import nnls, lsq_linear
sol_nnls = lsq_linear(eq.A, -eq.r, bounds=(0, np.inf), tol=1e-10).x
plt.figure(figsize=(4, 4))
plt.plot(sol_nnls, y[-1], '.')
plt.xlabel("Constrained inear regression solution")
plt.ylabel("Steady-state of Lotka Volterra model")
from scipy.stats import spearmanr
corr = spearmanr(sol_nnls, y[-1])
plt.title(f"Correlation: {corr.correlation:.2f}")
No description has been provided for this image
Out[162]:
Text(0.5, 1.0, 'Correlation: 0.81')
No description has been provided for this image
In [ ]:
 

What sets the hardness of linear regression?¶

Given linear regression with an interaction matrix $A$ and a vector of growth rates $r$, we can write the update rule as

In [ ]:
 
In [ ]:
 

Functional reduncancy is ill conditioning¶

  • Unlike $\epsilon$, which is a property of our particular model for generating interaction matrices $A$, the condition number is a property we can measure from any interaction $A$.

  • Two nearly-parallel rows of $A$ is equivalent to two species being nearly-redundant (interchangeable in the community).

In [ ]:
tmax = 1e8

n_val = 200
ic = np.random.uniform(size=n_val)
settling_times = []
eps_vals = np.logspace(-5, 0, 10)
condition_numbers = []
for eps in eps_vals:
    print(f"eps: {eps}")
    eq = RandomLotkaVolterra(n_val, random_state=0, eps=eps,  sigma=2.0, d=4.5, kfrac=0.2, connectivity=0.05)
    t, y = eq.integrate(tmax, ic)
    tind = compute_settling_time(y)
    settling_times.append(t[tind])
    condition_numbers.append(np.linalg.cond(eq.A))


plt.figure(figsize=(5, 5))
plt.loglog(eps_vals, settling_times, '.')
plt.xlabel("Distance to Exact Functional Redunancy")
plt.ylabel("Settling Time")

plt.figure(figsize=(5, 5))
plt.loglog(eps_vals, condition_numbers, '.k')
plt.xlabel("Distance to Exact Functional Redunancy")
plt.ylabel("Condition Number")
Out[ ]:
Text(0, 0.5, 'Condition Number')
No description has been provided for this image
No description has been provided for this image
In [7]:
tmax = 1e8

n_val = 200
ic = np.random.uniform(size=n_val)
settling_times = []
eps_vals = np.logspace(-5, 0, 10)
condition_numbers = []
for i in range(10):
    progress_bar(i, 20)
    eq = RandomLotkaVolterra(n_species=n_val, random_state=i, eps=np.random.exponential(scale=1e-3), sigma=2.0, d=4.5, kfrac=0.2, connectivity=0.05)
    # eq = GaussianLotkaVolterra(n_val, random_state=seed)
    # eq.A = 1.5 * normal_generalized(size=(n_val, n_val), beta=np.random.exponential(scale=0.5)) - d * np.identity(n_val)
    # eq.r = normal_generalized(size=n_val)
    t, y = eq.integrate(tmax, ic)
    tind = compute_settling_time(y)
    settling_times.append(t[tind])
    condition_numbers.append(np.linalg.cond(eq.A))

plt.figure(figsize=(5, 5))
plt.loglog(condition_numbers, settling_times, '.')
plt.xlabel("Condition Number")
plt.ylabel("Settling Time")
[#########           ] 
Out[7]:
Text(0, 0.5, 'Settling Time')
No description has been provided for this image

Krylov bound on linear solvers¶

The Krylov bound on linear solvers expresses the convergence time in terms of the condition number of the matrix.

$$ \tau = $$

In [10]:
## Loads a figure from the resources folder
Image("./resources/fig_overview.jpg", width=1200)
Out[10]:
No description has been provided for this image
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 

How does the settling time depend on the number of species?¶

In [76]:
n_vals = np.unique(np.logspace(1, 4, 200).astype(int)) # Range of species numbers to consider
tmax = 20000 # Maximum integration time
d = 21.5 # Density-limitation


## Loop over species numbers and compute the settling time for each
all_t = []
all_y = []
settling_times = []
for i,n in enumerate(n_vals):
    progress_bar(i, len(n_vals))

    eq = GaussianLotkaVolterra(n)
    eq.A = np.random.normal(size=(n, n)) - d * np.identity(n)
    eq.r = np.random.normal(size=n)
    ic = np.random.uniform(size=n)
    t, y = eq.integrate(tmax, ic)
    all_t.append(t.copy())
    all_y.append(y.copy())
    
    tind = compute_settling_time(y)
    if not np.isnan(tind):
        settling_times.append(t[tind])
    else:
        settling_times.append(np.nan)
[################### ] 

In [79]:
plt.loglog(n_vals, settling_times, '.k')
plt.xlabel("Number of species")
plt.ylabel("Settling time")
plt.xlim(None, 1000)
plt.ylim(1e1, None)
Out[79]:
(10.0, np.float64(1602.7528342389994))
No description has been provided for this image

Supertransients arise in hard optimization problems¶

  • Chimeras in the Kuramoto model
  • Pipe turbulence
  • Couple map lattices
  • Continuous time KSAT solvers
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 

How do transients depend on the initial conditions?¶

In [891]:
U, S, V = np.linalg.svd(eq.A)
# P = U[:, :10].T
# P1 = U[:, -10:].T
P_fast, P_slow = V[:10, :], V[-10:, :]
# P_slow = P_slow.T
def instantaneous_fli(y, t, dt, jac):
    
    # P2 = P_slow @ jac(t, y) @ P_slow.T

    # # ## Finite difference approximation of Jacobian
    eps = 1e-8
    jac_fd = np.zeros((P_slow.shape[0], P_slow.shape[0]))
    for i in range(jac_fd.shape[0]):
        pt0, pt1 = y.copy(), y.copy()
        pt1[i] += eps
        jac_fd[:, i] = (P_slow @ (eq(0, pt1) - eq(0, pt0))) / (P_slow @ (pt1 - pt0))
        # jac_fd[:, i] = (P_fast @ (eq(0, pt1) - eq(0, pt0))) / (P_fast @ (pt1 - pt0))
    jac_fd = np.array(jac_fd)
    P2 = jac_fd

    # Finite‐difference reduction:
    
    # eps = 1e-8
    # n_slow = P_slow.shape[1]   # should be 10
    # J_reduced_fd = np.zeros((n_slow, n_slow))
    # f0 = eq(0, y)  # baseline 100‐vector
    # for j in range(n_slow):
    #     p_j = P_slow[:, j]              # 100‐vector (slow direction j)
    #     y_perturbed = y + eps * p_j     # perturb along the jth slow basis
    #     f1 = eq(0, y_perturbed)         # 100‐vector
    #     # directional derivative of f in direction p_j:
    #     df_dp = (f1 - f0) / eps         # 100‐vector ≈ J_full @ p_j
    #     # project onto slow‐coordinates on the left:
    #     J_reduced_fd[:, j] = P_slow.T @ df_dp   # 10‐vector
    # P2 = J_reduced_fd

    # Let U[k, :] be the k-th slow‐mode represented as a 100‐vector.
    # U = P_slow.T       # shape (10, 100)
    # U = P_slow.T
    # J_full = jac(t, y)
    # # print(J_full.shape, U.shape)
    # # Preallocate the “scaled‐projection” matrix:
    # #   jac_fd[k, j] = (U[k, :] · (J_full[:, j])) / U[k, j]
    # n_modes, dim = U.shape
    # jac_fd  = np.zeros((n_modes, dim))

    # numerator = U @ J_full  # (10, 100)
    # # We only need the first n_modes columns (j=0..n_modes-1) to match the original FD over 10 coordinates.
    # # Denominator = U[:, :n_modes] (shape: 10 × 10)
    # num_slice = numerator[:, :n_modes]  # (10, 10)
    # den_slice = U[:, :n_modes]       # (10, 10)
    # # Compute jac_fd[k, j] = (u_k^T · (J_full · e_j)) / u_{k,j} for j=0..n_modes-1
    # jac_fd = num_slice / den_slice  # shape: (10, 10)
    # P2 = jac_fd

    

    Q = np.linalg.inv(np.eye(P2.shape[0]) - dt * P2)
    eigenvalues = np.linalg.eigvalsh(Q.T @ Q)
    fli_val = np.sqrt(np.max(eigenvalues))
    return fli_val


def find_fli(sol, tvals, jac):
    dtvals = np.diff(tvals)
    dtvals = np.hstack([dtvals, dtvals[-1]])
    fli = -np.inf
    for yval, dt, tval in zip(sol, dtvals, tvals):
        fli = max(fli, instantaneous_fli(yval, tval, dt, jac))
    return fli

all_fli = []
for sol, tvals in all_sol:
    all_fli.append(find_fli(sol, tvals, eq.jac))

plt.plot(n1_val, all_fli);
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[891], line 74
     72 all_fli = []
     73 for sol, tvals in all_sol:
---> 74     all_fli.append(find_fli(sol, tvals, eq.jac))
     76 plt.plot(n1_val, all_fli);

Cell In[891], line 69, in find_fli(sol, tvals, jac)
     67 fli = -np.inf
     68 for yval, dt, tval in zip(sol, dtvals, tvals):
---> 69     fli = max(fli, instantaneous_fli(yval, tval, dt, jac))
     70 return fli

Cell In[891], line 16, in instantaneous_fli(y, t, dt, jac)
     14     pt0, pt1 = y.copy(), y.copy()
     15     pt1[i] += eps
---> 16     jac_fd[:, i] = (P_slow @ (eq(0, pt1) - eq(0, pt0))) / (P_slow @ (pt1 - pt0))
     17     # jac_fd[:, i] = (P_fast @ (eq(0, pt1) - eq(0, pt0))) / (P_fast @ (pt1 - pt0))
     18 jac_fd = np.array(jac_fd)

KeyboardInterrupt: 
In [892]:
# n_res = 120
# n_res = 100
# n_res = 10
n_res = 10
# n_res = 10
# n_res = 25
# n_res = 10
# n_res = 5
# x = np.linspace(0.2, 0.202, n_res)
# y = np.linspace(0.2, 0.202, n_res)
x = np.linspace(0.2, 0.3, n_res)
y = np.linspace(0.2, 0.3, n_res)
x = np.linspace(0.0, 1.0, n_res)
y = np.linspace(0.0, 1.0, n_res)
xx, yy = np.meshgrid(x, y)
xx, yy = xx.ravel(), yy.ravel()
pts = np.vstack([xx, yy]).T

## Fixed disorder in initial conditions
np.random.seed(0)
# ic0 = np.random.uniform(size=eq.n)
icax = np.random.uniform(size=eq.n) * np.max(sol_final)
icbx = np.random.uniform(size=eq.n) * np.max(sol_final)
icay = np.random.uniform(size=eq.n) * np.max(sol_final)
icby = np.random.uniform(size=eq.n) * np.max(sol_final)



all_fli = []
for i, ic_val in enumerate(pts):
    progress_bar(i, len(pts))
    ic = bilinear_interpolation(icax, icbx, icay, icby, ic_val[0], ic_val[1])
    ic[ic < 0] = 0.0
    # print(ic[0])
    t, y = eq.integrate(tmax, ic, tolerance=1e-10)
    fli = find_fli(y, t, eq.jac)
    all_fli.append(fli)
    
[################### ] 

In [889]:
ms = 1200 * (20 / n_res)**2
plt.figure(figsize=(8,8))
all_fli = np.array(all_fli).astype(float)
vmin = np.percentile(all_fli, 5)
vmax = np.percentile(all_fli, 95)
plt.scatter(pts[:, 0], pts[:, 1], c=all_fli, s=ms, vmin=vmin, vmax=vmax)
# plt.colorbar()
plt.axis("off")
Out[889]:
(np.float64(-0.05), np.float64(1.05), np.float64(-0.05), np.float64(1.05))
No description has been provided for this image
In [867]:
 
Out[867]:
(np.float64(-0.05), np.float64(1.05), np.float64(-0.05), np.float64(1.05))
No description has been provided for this image
In [ ]:
    for pt, dt, tval in zip(fsol.y.T, dtvals, fsol.t):
        
        ## Finite difference approximation of Jacobian
        eps = 1e-8
        jac_fd = np.zeros((P_slow.shape[0], P_slow.shape[0]))
        for i in range(jac_fd.shape[0]):
            pt0 = pt.copy()
            pt1 = pt.copy()
            pt1[i] += eps
            # jac_fd[:, i] = (P_slow @ (eq(0, pt1) - eq(0, pt0))) / (P_slow @ (pt1 - pt0))
            jac_fd[:, i] = (P_fast @ (eq(0, pt1) - eq(0, pt0))) / (P_fast @ (pt1 - pt0))
        jac_fd = np.array(jac_fd)
        P2 = jac_fd

        # P2 = eq.jac(tval, pt)

        Q = np.linalg.inv(np.eye(P2.shape[0]) - dt * P2)
        eigenvalues = np.linalg.eigvalsh(Q.T @ Q)
        fli_val = np.sqrt(np.max(eigenvalues))

        fli = max(fli, fli_val)

    # project onto slow subspace
    sol = P_fast @ fsol.y
In [11]:
## Loads a figure from the resources folder
Image("./resources/fig_chaos.jpg", width=1200)
Out[11]:
No description has been provided for this image
In [ ]:
 
In [ ]:
 

Dimensionality reduction as preconditioning¶

  • What if we project the dynamics into a lower-dimensional space?

  • We start from our dynamics in the full $N$-dimensional space

$$ \dot{\mathbf{N}} = \mathbf{f}(\mathbf{N}) $$

  • We can instead study the projected by finding a projection operator $P \in \mathbb{R}^{N \times k}$, $k \ll N$, that projects the dynamics into a lower-dimensional space of dimension $k$

    $$ \mathbf{x} \equiv P \mathbf{N}, \quad \dot{\mathbf{x}} = P \dot{\mathbf{N}} $$

    where $\mathbf{x}(t) \in \mathbb{R}^k$ is a lower-dimensional state vector.

In [17]:
## Start by simulating a full-dimensional state vector

n_val = 200
eq = RandomLotkaVolterra(100, random_state=0)
ic = np.random.uniform(size=eq.n)
t, y = eq.integrate(tmax, ic)

plt.figure()
plt.semilogx(t, y);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, None)
Out[17]:
(0.01, np.float64(794328234.7242821))
No description has been provided for this image

How to choose a good projection $P$?¶

  • Since ill-conditioning leads to timescale separation, we seek a projection $P$ satisfying the following properties:

$$ \kappa(P A) \ll \kappa(A) $$

  • In numerical analysis, this is known as a preconditioner. If we have a stiff or ill-conditioned linear system, we can precondition it with an appropriate choice of $P$.

  • In PDEs, spectral methods can be thought of as a preconditioner. Some problems are ill-conditioned in the physical space, well-conditioned in the Fourier space.

In [64]:
print(f"The condition number of the full system is {np.linalg.cond(eq.A)}")
The condition number of the full system is 19686952.699464306
In [ ]:
 
In [ ]:
 

Singular value decomposition¶

  • Every matrix $A$ can be decomposed into a product of three matrices:

$$ A = U \Sigma V^\top $$

  • where $U \in \mathbb{R}^{N \times N}$ and $V \in \mathbb{R}^{N \times N}$ are orthogonal matrices (rotation matrices), and $\Sigma \in \mathbb{R}^{N \times N}$ is a diagonal matrix of positive singular values (in decreasing order).

  • We view $U$ as a rotation matrix that rotates our problem into more favorable coordinates.

In [69]:
U, S, V = np.linalg.svd(eq.A)

plt.figure()
plt.plot(S);
plt.xlabel("Rank")
plt.ylabel("Singular value")
plt.title("Singular values of the interaction matrix")
plt.semilogy();
No description has been provided for this image

How do we interpret the orthogonal matrices $U$ and $V$?¶

  • The right singular vectors $V_i$ isolate modes with dynamical timescales $\tau_i \sim 1 / \sigma_i$.
In [22]:
## Loads a figure from the resources folder
Image("./resources/fig_svd.jpg", width=1200)
Out[22]:
No description has been provided for this image
In [70]:
U, S, V = np.linalg.svd(eq.A)
P_fast, P_slow = V[:10, :], V[-10:, :]

yp_fast = P_fast @ y.T # First k modes
yp_slow = P_slow @ y.T # Last k modes


plt.figure()
plt.semilogx(t, y);
plt.xlabel("Time")
plt.ylabel("Species amplitude")
plt.xlim(1e-2, None)


plt.figure()
plt.semilogx(t, yp_fast.T);
plt.xlabel("Time")
plt.ylabel("Fast mode amplitude")
plt.xlim(1e-2, None)

plt.figure()
plt.semilogx(t, yp_slow.T);
plt.xlabel("Time")
plt.ylabel("Slow mode amplitude")
plt.xlim(1e-2, None)
plt.xlabel("Time")
plt.ylabel("Slow mode amplitude")
plt.xlim(1e-2, None)
Out[70]:
(0.01, np.float64(2542579.412082511))
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Preconditioning: We can compute the condition number of the dynamics projected into the subspaces spanned by the fast and slow modes.

$$ A_\text{slow} = P_\text{slow}^\top A P_\text{slow} $$

For a linear dynamical system, this would also correspond to the Jacobian of the projected system.

In [26]:
print("Original condition number: ", np.linalg.cond(eq.A))

print("Fast modes condition number: ", np.linalg.cond(P_fast @ eq.A @ P_fast.T))
print("Slow modes condition number: ", np.linalg.cond(P_slow @ eq.A @ P_slow.T))
Original condition number:  4429801842.761199
Fast modes condition number:  1.5851693011189856
Slow modes condition number:  153.51209647795702

Why does this work? The condition number alternatively can be expressed as the ratio of the largest and smallest singular values of the matrix.

$$ \kappa(A) = \frac{\sigma_{\text{max}}(A)}{\sigma_{\text{min}}(A)} $$

How general is this? Check out Weyl's inequalities discussed in Thibeault et al. 2024.









When would ill-conditioning emerge in a random ecosystem?¶

  • We revisit the case where the dynamics are not ill-conditioned, the Gaussian Lotka-Volterra model.
In [150]:
eq = GaussianLotkaVolterra(
    n=100, 
    d=8.5,
    random_state=0,
)
tmax = 1000
ic = np.random.uniform(size=eq.n)

t, y = eq.integrate(tmax, ic)

plt.figure()
plt.semilogx(t, y);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, t[-1])


print(f"The number of surviving species is {np.sum(y[-1] > 1e-2)} out of {eq.n}")
The number of surviving species is 49 out of 100
No description has been provided for this image
  • We now perform a simple evolutionary algorithm to evolve the ecosystem to have more surviving species.

  • We start with a random interaction matrix, simulate the ecosystem, and find the surviving species.

  • If a species is excluded from the ecosystem, we replace it in the next generation with a new species with interactions averaged from two randomly chosen surviving species.

  • If $i$ is the index of the extinct species, $j$ and $k$ are the indices of the two surviving species, and $a_{ij}$ is the interaction between species $i$ and $j$, we update the interaction matrix as

    $$ a_{ik} = 0.05 a_{ij} + 0.95 a_{ik} $$

  • We repeat this process for 100 generations.

In [151]:
all_fitnesses = []
all_condition_numbers = []
for j in range(100):

    ## Simulate the ecosystem with the current interaction matrix
    t, y = eq.integrate(tmax, ic)
    sol = y[-1]
    
    ## Find the surviving species
    where_survive = (sol > 1 / eq.n)
    surviving_inds = np.where(where_survive)[0] 
    extinct_inds = np.setdiff1d(np.arange(eq.n), surviving_inds)


    ## Recombination: replace extinct species with new species averaged from surviving species
    new_matrix = eq.A.copy()
    for i in extinct_inds:
        parent_ind1, parent_ind2 = np.random.choice(surviving_inds, size=2, replace=False)
        interaction_ind = np.random.choice(np.setdiff1d(np.arange(eq.n), [parent_ind1, parent_ind2, i]))
        new_matrix[i] = eq.A[i]
        new_matrix[i, interaction_ind] = 0.5 * eq.A[parent_ind1, interaction_ind] + 0.5 * eq.A[parent_ind2, interaction_ind]

    ## Update the interaction matrix
    eq.A = new_matrix.copy()
    fitness = np.sum(where_survive)
    print(f"Iteration {j}: The number of surviving species is {fitness} out of {eq.n}")

    ## Save the fitness and condition number for this generation
    all_fitnesses.append(np.sum(where_survive))
    all_condition_numbers.append(np.linalg.cond(eq.A[surviving_inds][:, surviving_inds]))
Iteration 0: The number of surviving species is 49 out of 100
Iteration 1: The number of surviving species is 50 out of 100
Iteration 2: The number of surviving species is 52 out of 100
Iteration 3: The number of surviving species is 53 out of 100
Iteration 4: The number of surviving species is 53 out of 100
Iteration 5: The number of surviving species is 54 out of 100
Iteration 6: The number of surviving species is 55 out of 100
Iteration 7: The number of surviving species is 55 out of 100
Iteration 8: The number of surviving species is 59 out of 100
Iteration 9: The number of surviving species is 59 out of 100
Iteration 10: The number of surviving species is 57 out of 100
Iteration 11: The number of surviving species is 56 out of 100
Iteration 12: The number of surviving species is 58 out of 100
Iteration 13: The number of surviving species is 59 out of 100
Iteration 14: The number of surviving species is 57 out of 100
Iteration 15: The number of surviving species is 59 out of 100
Iteration 16: The number of surviving species is 60 out of 100
Iteration 17: The number of surviving species is 60 out of 100
Iteration 18: The number of surviving species is 60 out of 100
Iteration 19: The number of surviving species is 60 out of 100
Iteration 20: The number of surviving species is 60 out of 100
Iteration 21: The number of surviving species is 60 out of 100
Iteration 22: The number of surviving species is 61 out of 100
Iteration 23: The number of surviving species is 61 out of 100
Iteration 24: The number of surviving species is 61 out of 100
Iteration 25: The number of surviving species is 61 out of 100
Iteration 26: The number of surviving species is 62 out of 100
Iteration 27: The number of surviving species is 62 out of 100
Iteration 28: The number of surviving species is 62 out of 100
Iteration 29: The number of surviving species is 62 out of 100
Iteration 30: The number of surviving species is 62 out of 100
Iteration 31: The number of surviving species is 62 out of 100
Iteration 32: The number of surviving species is 58 out of 100
Iteration 33: The number of surviving species is 58 out of 100
Iteration 34: The number of surviving species is 58 out of 100
Iteration 35: The number of surviving species is 63 out of 100
Iteration 36: The number of surviving species is 64 out of 100
Iteration 37: The number of surviving species is 68 out of 100
Iteration 38: The number of surviving species is 68 out of 100
Iteration 39: The number of surviving species is 68 out of 100
Iteration 40: The number of surviving species is 68 out of 100
Iteration 41: The number of surviving species is 69 out of 100
Iteration 42: The number of surviving species is 70 out of 100
Iteration 43: The number of surviving species is 70 out of 100
Iteration 44: The number of surviving species is 70 out of 100
Iteration 45: The number of surviving species is 70 out of 100
Iteration 46: The number of surviving species is 70 out of 100
Iteration 47: The number of surviving species is 58 out of 100
Iteration 48: The number of surviving species is 64 out of 100
Iteration 49: The number of surviving species is 63 out of 100
Iteration 50: The number of surviving species is 69 out of 100
Iteration 51: The number of surviving species is 69 out of 100
Iteration 52: The number of surviving species is 67 out of 100
Iteration 53: The number of surviving species is 66 out of 100
Iteration 54: The number of surviving species is 67 out of 100
Iteration 55: The number of surviving species is 70 out of 100
Iteration 56: The number of surviving species is 70 out of 100
Iteration 57: The number of surviving species is 70 out of 100
Iteration 58: The number of surviving species is 70 out of 100
Iteration 59: The number of surviving species is 71 out of 100
Iteration 60: The number of surviving species is 71 out of 100
Iteration 61: The number of surviving species is 71 out of 100
Iteration 62: The number of surviving species is 72 out of 100
Iteration 63: The number of surviving species is 72 out of 100
Iteration 64: The number of surviving species is 72 out of 100
Iteration 65: The number of surviving species is 72 out of 100
Iteration 66: The number of surviving species is 73 out of 100
Iteration 67: The number of surviving species is 74 out of 100
Iteration 68: The number of surviving species is 75 out of 100
Iteration 69: The number of surviving species is 75 out of 100
Iteration 70: The number of surviving species is 75 out of 100
Iteration 71: The number of surviving species is 76 out of 100
Iteration 72: The number of surviving species is 77 out of 100
Iteration 73: The number of surviving species is 78 out of 100
Iteration 74: The number of surviving species is 79 out of 100
Iteration 75: The number of surviving species is 80 out of 100
Iteration 76: The number of surviving species is 80 out of 100
Iteration 77: The number of surviving species is 81 out of 100
Iteration 78: The number of surviving species is 81 out of 100
Iteration 79: The number of surviving species is 81 out of 100
Iteration 80: The number of surviving species is 81 out of 100
Iteration 81: The number of surviving species is 80 out of 100
Iteration 82: The number of surviving species is 80 out of 100
Iteration 83: The number of surviving species is 81 out of 100
Iteration 84: The number of surviving species is 81 out of 100
Iteration 85: The number of surviving species is 81 out of 100
Iteration 86: The number of surviving species is 81 out of 100
Iteration 87: The number of surviving species is 81 out of 100
Iteration 88: The number of surviving species is 81 out of 100
Iteration 89: The number of surviving species is 82 out of 100
Iteration 90: The number of surviving species is 82 out of 100
Iteration 91: The number of surviving species is 84 out of 100
Iteration 92: The number of surviving species is 85 out of 100
Iteration 93: The number of surviving species is 85 out of 100
Iteration 94: The number of surviving species is 85 out of 100
Iteration 95: The number of surviving species is 85 out of 100
Iteration 96: The number of surviving species is 85 out of 100
Iteration 97: The number of surviving species is 85 out of 100
Iteration 98: The number of surviving species is 85 out of 100
Iteration 99: The number of surviving species is 85 out of 100
In [152]:
plt.plot(all_fitnesses);
No description has been provided for this image
In [155]:
plt.semilogy(all_condition_numbers)
plt.xlabel("Generation")
plt.ylabel("Condition number")
Out[155]:
Text(0, 0.5, 'Condition number')
No description has been provided for this image
In [12]:
## Loads a figure from the resources folder
Image("./resources/fig_evolve.jpg", width=1200)
Out[12]:
No description has been provided for this image
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
tmax = 1e8 ## Maximum integration time

eq = RandomLotkaVolterra(
    n_species=n_val, # Number of species
    sigma=2.0, # Amplitude of interactions
    d=4.5, # Density-limitation strength
    kfrac=0.5, # Fraction redundant interactions
    eps=0.00005, # Functional redundancy
    random_state=0, # Random seed
    connectivity=0.05 # Connectivity of network
)

ic = 2*np.random.uniform(size=eq.n)
tlim = (0, 50000)
t, y = eq.integrate(tmax, ic)
plt.figure()
plt.semilogx(t, y, color="r", linewidth=0.8);

ic += 1e-3
t, y = eq.integrate(tmax, ic)
plt.semilogx(t, y, color="k", linewidth=0.8);

plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, t[-1])

# settling_time = compute_settling_time(y)
# print(f"Settling time: {t[settling_time]:.2f}")
In [ ]:
n1_val = np.linspace(0, 1, 100)
n1_val = np.linspace(0.2, 0.3, 100)
n1_val = np.linspace(0.2, 0.21, 100)
icp = ic.copy()
all_sol = []
for i, n1 in enumerate(n1_val):
    progress_bar(i, len(n1_val))
    icp[0] = n1
    t, y = eq.integrate(tmax, icp)
    all_sol.append([y, t])



    
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 

The ill-conditioned Lotka-Volterra model¶

We consider random ecosystems given by the generalized Lotka-Volterra equation,

$$ \frac{dN_i}{dt} = N_i \left( r_i + \sum_{j=1}^N A_{ij} N_j \right) $$

where $N_i$ is the population of species $i$, $r_i$ is the intrinsic growth rate of species $i$, and $A_{ij}$ is the interaction coefficient between species $i$ and $j$. The steady-state solutions of this equation has the form

$$ -A \mathbf{N}^* = \mathbf{r} $$

where $N_i \geq 0$ for all $i$.

In this notebook, we explore the behavior of this model when the interaction matrix $A_{ij}$, and the growth rates $r_i$ are drawn from random distributions. Specifically, we consider the case where $r_i \sim \mathcal{N}(0,1)$, and the matrix $A$ has the form

$$ A = P^\top (Q - d\, I) P + \epsilon E $$

where $Q_{ij} \sim \mathcal{N}(0,1)$, $E_{ij} \sim \mathcal{N}(0,1)$, $P$ is a low-rank matrix imposing functional redundancy, $d$ is a constant density-limitation, and $\epsilon \ll 1$ is a small constant.

In [8]:
from scipy.integrate import solve_ivp

## Load local functions
import sys
import base
import utils
from base import RandomLotkaVolterra
from utils import levenshtein


tlim = (0, 20000) # Integration interval
n_val = 200 # Number of species

## Initialize the model
eq = RandomLotkaVolterra(n_val, sigma=2.0, d=4.5, kfrac=0.2, eps=0.001, random_state=0, connectivity=0.01)

## Initial conditions
ic = np.random.uniform(size=eq.n)

## Numerical integration
fsol = solve_ivp(eq, tlim, ic[:n_val], jac=eq.jac, **eq.integrator_args)


## Check that the found solution is stable
sol_final = fsol.y[:, -1]
jac_final = eq.jac(0, sol_final)
eigs = np.linalg.eigvals(jac_final)
print(f"Numerical stability observed: {np.all(np.real(eigs) < 0)}")

## Plot the solution
plt.figure(figsize=(7, 4))
plt.semilogx(fsol.t, fsol.y.T, color="k", lw=1.5, alpha=0.3);
plt.xlabel("Time")
plt.ylabel("Species abundance")
plt.xlim(1e-2, np.max(fsol.t))
plt.ylim(0, None)

eigvals = np.linalg.eigvals(jac_final)
plt.figure(figsize=(4, 4))
plt.plot(np.real(eigvals), np.imag(eigvals), '.k')
plt.xlabel("Real part")
plt.ylabel("Imaginary part")
plt.axis("equal")
plt.title("Eigenvalues of Jacobian")
Numerical stability observed: True
Out[8]:
Text(0.5, 1.0, 'Eigenvalues of Jacobian')
No description has been provided for this image
No description has been provided for this image
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: